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(I.)  Introduction 


Grant  No.  AFOSR-81 -01 70  has  an  objective  which  is  well- 
summarized  by  the  Grant  title:  "Feasibility  studies  of  optical 
processing  for  image  bandwidth  compression  schemes."  It  is  the 
intent  of  research  sponsored  under  this  Grant  to  direct  investi¬ 
gation  into  the  following  issues: 

(a)  formulation  of  alternative  archi techtural  concepts 
for  image  bandwidth  compression,  i.e.,  the  formula¬ 
tion  of  components  and  schematic  diagrams  which 
differ  from  conventional  di qi tal  bandwidth  compres¬ 
sion  schemes  by  being  implemented  by  various 
optical  computation  methods; 

(b)  simulation  of  optical  processing  concepts  for  image 
bandwidth  compression,  so  as  to  gain  insight  into 
typical  performance  parameters  and  elements  of  system 
performance  sensitivity; 

(c)  maturation  of  optical  processing  for  image  band¬ 
width  compression  until  the  overall  state  of  optical 
methods  in  image  compression  becomes  equal  to  that 
of  digital  image  compression. 

It  is  the  last  of  these,  item  (c),  which  represents  the 
continuing  strategic  objective  of  the  efforts  being  carried  on 
under  Grant  No.  AF0SR-81-01 70.  It  is  important  to  remember  that 
the  major  attention  given  to  image  bandwidth  compression  has 
been  for  methods  most  conveniently  implemented  by  digital  compu- 
tions.  As  flexible  and  multipurpose  are  digital  methods,  there 


may  always  be  operational  circumstances,  environments,  or  con¬ 
straints  where  the  availability  of  a  different  technology  is 
important.  However,  with  the  concentration  upon  digital  compu¬ 
tations,  which  has  characterized  most  research  on  bandwidth 
compression,  alternative  methods  in  optics  have  suffered.  Thus 
the  purpose  of  research  sponsored  under  this  Grant  is  to  serve 
as  a  source  of  alternatives  for  future  concepts  in  bandwidth 
compression,  so  that  the  environment  for  compression  technology 
need  not  be  dominated  by  one  methodology. 
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(II.)  Overview  of  the  Report 


The  research  currently  sponsored  under  Grant  No.  AFOSR-81- 
0170  consists  of  several  distinct  and  separate  activities.  The 
separate  research  efforts  are  unified  by  a  common  theme:  the 
application  of  optical  processing  for  image  bandwidth  compres¬ 
sion.  Within  this  common  theme,  however,  the  separate  research 
projects  are  not  completely  related  to  each  other.  There¬ 
fore,  this  report  is  put  together,  literally,  as  a  number  of 
independent  reports.  The  separate  sections  of  the  report,  which 
follow  this  section,  are  intended  to  be  read  separately  and  in¬ 
dependently  of  any  other  section.  Each  section  has  its  own 
references  and  its  own  figure  labellings,  for  example. 

The  separate  sections  of  this  report,  and  the  research 
problems  dealt  with  in  each  section,  are  summarized  in  the  fol¬ 
lowing" 


(1)  Data  compression  by  optical  tomography,  with  data 
reconstruction  by  optical  convolution  and  back  pro¬ 
jection  (see  Section  III). 

(2)  Adaptive  data  compression  by  modification  of  a  previously 
demonstrated  technique,  IDPCM,  to  an  efficient  spatially 
recursive  scheme  (see  Section  IV). 


(3)  Adaptive  image  processing  by  using  spatial  transformations 
to  create  a  spatially  stationary  image  (see  Section  V). 


Tomography  and  the  Projection  Matrix 


Tomography  is  a  procedure  which  decomposes  a  two-dimensional  image 
into  a  series  of  one-dimensional  projections,  each  made  at  a  different 
angle  through  the  original  image.  A  projection  is  produced  by  integrating 
the  image  data  in  one  direction  across  the  image.  Along  the  horizontal 
axis,  for  example  a  projection  is  defined  by: 

oo 

Pe(x)  *  j  f(x,y)dy 

—  6* 

where  P0  is  the  projection  at  angle  ®  and  f(x,y)  is  the  original  image. 
Subsequently,  either  the  coordinate  system  or  the  object  is  rotated  and  the 
next  projection  is  calculated. 

In  digital  tomography,  the  summed  (integrated)  data  results  are 
registered  by  a  or.e- dimensional  array  of  discrete  sensors  or  detectors. 

This  string  of  detectors  must  be  large  enough  to  record  all  of  the  data  at 
each  of  the  possible  angles.  For  a  square  image,  the  maximum  number  of 
detectors  is  required  at  the  angles  of  45°  and  135°;  see  Figure  1.  At 
other  angles  however,  the  ends  of  the  detector  array  swing  outside  of  the 
image,  thus  registering  artificial  zero  data.  The  scalloped  ends  of  the 
projection  matrix,  visible  in  Figure  2,  are  the  result  of  these  artificial 
zeros. 

The  projection  matrix  of  any  image  displays  certain  consistent  char¬ 
acteristics.  As  seen  in  the  projection  matrix  of  Figure  2,  there  is  a 
sinusoidal  design  woven  into  the  matrix;  this  pattern  is  present  in  the 
projection  matrix  derived  from  any  image. 

Each  horizontal  line  in  the  projection  matrix  corresponds  to  the 
projection  data  gathered  by  all  N  detectors  at  one  specific  angle;  a 
horizontal  line,  therefore,  is  referred  to  as  a  data  sequence  from  "within" 


Figure  1:  demonstrates  how  the  number  of  detectors  necessary  varies  as  a 
function  of  projection  angle.  At  a  projection  angle  of  45°  (a)  seven 
detectors  are  required,  but  at  a  projection  angle  of  0°  (b),  only  five 
detectors  are  reeded;  the  extra  detectors  at  the  endpoints  swing  out  of 
the  image. 


Figure  2:  Projection  matrix  or  sinogram 


.  •.  \  \  V  V  *,  *.  v  \  .  --  V  *.  *-  -  - 
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a  projection  (See  Figure  3).  Conversely,  a  vertical  line  through  the 
projection  matrix  represents  the  data  integrated  at  one  detector  as  it 
travels  through  all  M  angles;  this  is  referred  to  as  data  "across"  pro¬ 
jections  (See  Figure  4). 

To  reconstruct  the  original  image,  filtered  back-projection  is  done. 
Each  projection  first  is  convolved  with  a  high  pass,  rho  filter.  Following 
this,  the  filtered  projections  are  spread  back  across  the  image  plane  at 
their  original  angle  and  then  the  results  from  all  angles  are  summed. 


Redundancy  in  a  Projection  Matrix 


The  purpose  of  data  compression  techniques  is  to  remove,  or  at 
least  reduce,  the  redundant  or  correlated  information,  thereby  retaining 
only  orthogonal  or  decorrelated  data.  It  is  well  known  that  the  data  in 
most  images  if  highly  redundant.  From  this  knowledge,  one  might  surmise 
that  the  projection  matrix  also  contains  a  great  deal  of  redundant  infor 
nation.  To  test  this,  a  method  was  devised  to  measure  the  entropy  in  a 
projection  matrix,  both  within  and  across  projections.  The  information 
or  entrcpy,  H,  can  be  calculated  using: 


He  1o92  F  bits 

Hd  *  £  I’d:  ,092  p„  bits 

i  **  1  d; 


(within  projections) 
(across  projections) 


where  n  is  the  number  of  cells  in  the  data  histogram, ©is  a  specific  angle, 
d  is  a  specific  detector  and  pQ-  or  pd.  is  the  probability  of  data  being 
contained  in  the  1th  cell  of  the  histogram. 


Figure  3:  The  values  of  the  original  image  (a),  integrated  along  y1  are 
sensed  by  each  of  the  detectors  along  x‘;  for  example,  the  result  at  detector 
number  three  is  the  sum  of  values  in  the  shaded  region  of  (a).  For  the 
projection  angle  0  ,  the  results  from  each  detector  are  placed  into  the 
the  projection  matrix  (b)  along  a  horizontal  line. 


Figure  4:  A  vertical  line  in  the  projection  matrix  (b),e.g.,  along  detector 
number  three,  represents,  in  the  original  image  (a),  the  integrated  result 
from  a  single  detector  (in  this  example,  number  three)  as  it  swings  through 
all  possible  angles. 


The  redundancy  in  a  line  of  projection  data  then  can  be  calculated: 

Re  =  Hmax  '  H©  (within)  or  Rd  =  Hmax  -  (across) 

where  is  the  maximum  amount  of  information  possible. 

Method 

The  original  image,  the  source  of  the  projection  data,  was  128x128 
pixels  in  size.  For  every  projection  angle,  182  detectors  (each  detector 
being  1  pixel  in  width)  were  used  and  a  total  of  100  projections  were  made. 
For  ease  in  calculation,  the  data  in  the  matrix,  with  an  original  range  from 
0  to  34323,  were  scaled  to  a  range  from  0  to  255.  A  histogram  of  data 
values  was  computed  for  each  line  in  the  matrix.  The  histogram  was  divided 
into  16  cells,  each  with  a  width  of  16  intensity  levels.  The  maximum 
amount  of  information  possible  would  be  present  if  all  cells  in  the  histo¬ 
gram  were  equally  probable: 

"max  =  l0S2  16  =  4  bUs 

Results 

Figure  5  shows  the  redundancy  within  each  projection  as  a  function 
of  projection  angle.  The  amount  of  redundancy  appears  to  be  strongly 
dependent  upon  the  angle  of  projection.  The  redundancy  is  highest  near 
the  angles  of  0°,90°  and  180°  and  is  lowest  near  angles  of  45°  and  135°. 

The  maximum  obtained  redundancy  of  2.184  (54.6%  of  the  maximum)  occurs  at 
an  angle  of  0°and  the  smallest  value  of  0.0577  (1.4%  of  maximum)  is  at 


Greater  redundancy  at  certain  angles  primarily  appears  to  be  a  function 


of  how  the  tomographic  process  operates  on  a  square  image.  Figure  6  shows 
why  such  an  effect  may  occur. 

Redundancy  across  projections  is  shown  in  Figure  7.  Detectors  near 
the.  endpoints  of  the  projection  matrix,  those  that  swing  out  of  the  image 
near  angles  of  0°,90°  and  18C°  were  excluded  from  analysis  to  avoid  the 
introduction  of  artificial  zero  data.  The  values  vary,  for  the  most  part 
between  redundancies  of  1  to  2  bits  (25-50%)  for  each  detector. 

Compression  of  Projections  using  DPCM 

Differential  pulse  code  modulation  (DPCM)  coders  operate  on  the 
principle  of  quantizing  prediction  error  values  rather  than  actual  data 
values.  Based  on  the  recent  history  of  the  signal,  a  predictor  £(n) 
is  made  to  approximate  the  actual  x(n).  Then,  rather  than  quantizing 
and  storing  x(n),  the  difference: 

d(n)=x(n)-x(n) 

is  quantized  and  stored.  Subsequently,  the  decoder  attempts  to  recover 
the  original  signal  by  essentially  integrating  the  quantized  differences 
between  samples.  Figure  8  shows  a  block  diagram  of  the  basic  closed  loop 
DPCM  system. 

For  a  data  sequence  x(n)  with  a  normalized  autocovariance,  ^ 1 , 
the  variance  of  the  difference  signal,  which  is  the  input  to  the  quantizer, 
is  significantly  smaller  than  that  of  the  original  data.  Since  the  quant¬ 
ization  error  variance  is  directly  proportional  to  the  variance  of  the:  quant¬ 
izer  input,  it  is  possible  to  lower  the  quantization  bit  rate  to  a  specific 
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Figure  5  demonstrates  how  redundancy  varies  as  a  function  of  projection  angle  due 
to  the  tomographic  process  operating  on  a  square  image.  Assuming  a  completely 
redundant  original  image,  the  projection  data  at£=0°  is  totally  redundant  also; 
however,  the  projection  at  @=45°  is  only  about  50%  redundant. 


level  and  maintain  the  same  signal  to  noise  ratio  (SNR).  Further  reductions 
in  the  bit  rate,  past  the  new  level,  also  are  possible  but  at  the  expense 
of  the  SNR. 

The  closed  loop  or  feedback-around  quantizer  diagram  of  Figure  8 
ensures  that  the  quantization  errors  do  not  accumulate  into  subsequent 
samples,  i.e.,  the  errors  should  be  independent. 

There  are  a  variety  of  prediction  schemes  which  can  be  used.  The 
two  used  in  this  study  are: 

1.  Normalized  autocovariance  or  rho  method 

2.  Average  rate  of  change  or  slope  method. 

Prediction  Methods 

For  both  of  the  prediction  methods,  coding  of  the  projection  matrix  is 
carried  out  one  line  at  a  time.  The  first  value  in  each  line  is  carried 
through  the.  system  at  its  full  bit  rate  and  the  remaining  values  in  the  line 
are  quantized  at  the  reduced  rate. 

The  prediction  errors,  the  input  to  the  quantizer  in  a  DPCM  system, 
tend  to  be  distributed  as  a  Laplacian  (exponential)  probability  density 
(Gray,  1983).  Decision  and  reconstruction  levels  have  been  derived  to 
minimize  the  mean  square  error  in  quantizing  a  Laplacian.  These  levels, 
given  in  Pratt  (1978),  are  used  in  the  quantizer  of  Figure  8.  The  mean 
square  error  between  the  original  and  decoded  values  of  the  projection 
matrix  were  computed  for  each  line  and  subsequently  were  averaged  across 
lines. 

In  order  to  perform  DPCM  on  lines  of  constant  sample  size,  the  scalloped 
ends  of  the  projection  matrix  were  trimmed  before  encoding,  producing  a 
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rectangular  matrix.  The  data  values  which  were  removed  were  those 
collected  by  the  detectors  farthest  from  the  origin,  those  that  swing 
out  of  the  image  at  angles  near  0°,  90°  and  180°.  Tomographic  reconstruction 
without  these  data  produces  a  circular  image. 

Rho  Prediction.  The  Rho  method  involves  a  linear  predictor  of  the 

form: 

a. 


x(n)  =  ^h.x(n-j) 


and  the  quantizer  input  for  a  first  order  predictor  is: 


d(n)  =  x(n)  -  h^x(n-l ) 


e.«“C  i  4  Ifrri 


Taking  the  square  and  expected  value  of  the  equation  above,  it  can  be  shown 

c-*  .  u-ah,e.  .  V) 

To  minimize  the  mean  square  error,  take  the  partial  derivative  of  the  variance 
with  respect  to  h^  and  set  it  equal  to  zero. 


C ^  -° 

cM. 


y*(  =  -Per  iniH'M'Wm 

Estimation  of  x  by  this  first  order  predictor  assumes  that  you 
have  no  other  information  concerning  the  expected  shape  of  the  projection 
data.  The  slope  method,  however,  relies  on  the  fact  that  some  general 
assumptions  can  be  made  concerning  the  projection  matrix. 

Slope  Method.  As  mentioned  earlier,  all  projection  matrices  share 
the  sinogram  pattern.  Based  on  this,  it  should  be  possible  to  make  a  more 
accurate  estimate  of  x  based  not  only  on  recent  signal  history  but  also  on 
where  x  lies  in  the  matrix  and  on  whether  DPCM  is  being  carried  out  within 
or  across  projections. 


The  trimmed  projection  matrix  is  partitioned  into  a  set  of  blocks, 
see  Figure  9.  For  each  block,  an  estimate  is  made  of  the  average  rate 
of  change  between  x(n)  and  x(n-l),  either  across  or  within  projections. 
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This  slope  value  is  then  used  as  a  predictor  in  both  the  coding  and  decoding 
stages: 

x\n)  *  x(n-l)  +  slope(m,i) 

where  x(n-l)  is  contained  in  the  block  (m,i)  of  the  matrix.  After  their 
calculation,  these  slope  values  are  retained  for  the  decoding  phase, 
consequently  adding  to  the  final  bit  rate  for  the  coded  matrix. 

Results 

Table  1  gives  the  mean  square  error  for  each  prediction  method  and 
for  various  bit  rates.  Before  compression,  the  trimmed  projection  matrix 
contained  values  ranging  from  a  minimum  of  6776  to  a  maximum  of  34323  with 
a  mean  of  21670.  Therefore,  a  mean  square  error  (MSE)  of  100,000.,  for 
example,  represents  an  average  error  of  317  per  element  of  the  matrix. 

For  the  minimum  value  this  represents  an  average  of  +  4.7%  error,  while 
for  the  maximum  value,  it  is  an  average  error  of  +  0.92%. 


Table  1 


a)  Within  Projections 


Rho 

Prediction 

Slope  Prediction 

MSE 

bit  rate 

MSE 

bit  rate 

54,745 

4.09 

41,681 

4.20 

115,784 

3.10 

74,831 

3.21 

270,901 

2.11 

158,945 

2.22 

986,390 

1.12 

483,260 

1.23 

%  *  *  •  •  »  -  o  o  * 
■  v.- 


Figure  9:  Trimmed  projection  matrix  which  has  been  partitioned  into  a 
set  of  blocks.  Each  block,  fcr  example,  may  be  10x10,  meaning  10  angles 
long  and  10  detectors  wide. 


Table  1  cont. 


b)  Across  Projections 

Rho  Prediction  _  Slope  Prediction 

MSE  bit  rate  MSE  bit  rate 


20,220 

4.13 

16,176 

4.23 

80,437 

3.14 

54,698 

3.25 

307,745 

2.14 

228,544 

2.25 

2,540,996 

1.15 

1,233,708 

1.26 

Once  the  compressed  projections  have  been  decoded,  the  next  step  is 
to  reconstruct  the  original  image  by  back-projection.  The  following 
photographs  show  the  image  in  its  original  and  reconstructed  forms  after 
being  subject  to  various  amounts  of  data  compression  and  different 
prediction  schemes. 


Figure  10 


a)  original  image 


b)  reconstructed  image, 
no  data  compression, 

16  bits/element  in  matrix 


c)  coding  within  projections 
4.09  bits/element  in  matrix 
rho  prediction  method 


d)  coding  within  projections 
4.2  bits/element  in  matrix 
slope  prediction  method 


A  post  hoc  analysis  of  the  error  distribution,  the  input  to  the 
quantizer,  was  done  to  determine  whether  the  assumption  of  a  Laplacian 
distribution  was  justified.  Figure  11  shows  the  error  [d(n)]  distributions 
resulting  from  both  rho  and  slope  prediction  methods.  The  rho  distribution 
contains  obvious  deviations  from  an  exponential,  but  the  slope  distr bution, 
by  visual  analysis,  appears  to  match  ar  exponential  quite  well.  The  main 
difference  between  the  slope  distribution  and  a  theoretical  Laplacian 
[p(x)=  j  exp(-aixi)]  occurs  at  the  peak  of  the  distribution  and  in  the  tails 
The  mode  of  the  obtained  distribution  is  lower  than  expected,  while  the 
tails  are  somewhat  higher.  Analysis  of  the  obtained  slope  distribution 
using  Pearson's  Chi  Square  Goodress  of  Fit  Statistic,  however,  indicates 
that  it  cannot  be  accurately  modeled  as  a  Laplacian  (p  >  0.995).  Based  on 
this  finding,  the  question  arises  as  to  whether  a  better  way  to  quantize 
the  errors  is  possible.  Preliminary  results  obtained  in  answering  that 
question  suggest  that  better  quantization  can  be  achieved  when  the  decision 
and  reconstruction  levels  are  extended  further  towards  the  tails  of  the 
distribution. 

Figure  12  shows  the  results  of  extending  the  quantizer  levels  outward. 
The  quantization  for  this  image  was  performed  manually,  by  visual  analysis 
or  the  histogram.  This  reconstruction,  with  a  bit  rate  of  2.1  bits  per 
element  of  the  projection  matrix,  appears  to  be  superior  to  the  image  of 
Figure  lOg  which  is  a  slightly  higher  bit  rate  (2.22  bits);  in  fact  Figure 
12  seems  to  be  almost  comparable  to  the  quality  of  some  images  quantized 
to  approximately  3  bits.  Although  a  sharper  image  is  obtained  by  adjusting 
the  quantization  levels,  a  large  amount  of  noise  still  appears  in  the 
reconstructed  image. 


Figure  12:  Reconstruction  from  projections  quantized  by  moving  decision 
and  reconstruction  levels  out  toward  the  tails  of  the  error  distribution 
Coding  done  within  projections,  bit  rate  of  2,1  bits/element  in  matrix 
ard  slope  prediction  method  was  used. 


Conclusion 


Tomographic  projections,  originally  16  bit  data,  can  be  compressed 
to  a  bit  rate  of  approximately  3.2  bits  per  element  of  projection  matrix 
(a  compression  ratio  of  5:1). and  still  produce  a  recognizable  image; 
however,  the  image  at  this  rate  is  of  poor  quality  and  unsuitable  for 
most  applications. 

An  improvement  in  results  can  be  achieved  by  adjusting  the  decision 
and  reconstruction  levels  of  the  quantizer.  Another  possible  avenue  for 
improvement  would  be  to  employ  more  sophisticated  prediction  procedures. 
However,  a  fundamental  problem  involved  in  this  technique  seems  to  rest 
upon  the;  fact  that  any  errors  in  the  decoded  projection  matrix  are 
amplified  by  the  filtering  done  before  back-projection.  This  high  pass 
filter  is,  most  likely,  the  cause  of  the  noisy  appearance  of  the  recon¬ 
structed  images.  Since  this  filter  is  an  integral  part  of  back-projecting 
tomographic  projections,  it  would  appear  that,  while  further  research 
into  this  area  may  improve  the  results  demonstrated  here,  the  amount  of 
improvement  possible  may  not  be  significant  enough  to  warrant  the  additional 
investigation. 


IV.)  Background  Discussion 


Image  data  compression  methods  can  be  classified  in  two 
basically  different  categories.  One  category  is  processing  in 
the  spatial  domain.  Another  is  "transform  coding".  In  the  first 
category  are  those  methods  which  exploit  redundancy  in  the  spatial 
data.  Redundancy  is  a  characteristic  which  is  related  to  pre¬ 
dictability,  randomness  in  the  data.  For  example,  an  image  of 
constant  gray  levels  is  fully  predictable  once  the  gray  level  of 
the  first  pixel  is  known.  On  the  other  hand,  a  white-noise  random 
signal,  such  as  that  seen  on  a  TV  screen  when  no  program  is  being 
broadcast,  is  totally  unpredictable  and  every  pixel  has  to  be 
stored  to  produce  the  image. 

Interpolated  difference  pulse-code  modulation,  developed  by 
Hunt  in  1977  [2],  is  a  successful  method  of  data  compression  in 
spatial  domain.  Since  then,  several  variations  of  this  method 
have  been  developed.  Recently,  a  significant  improvement  to 
IDPCM,  named  recursive  IDPCM,  was  demonstrated  by  Hunt  and  Cao  [1]. 

Recursive  IDPCM  has  two  main  features: 

(1)  This  is  a  very  efficient  image  data  compression  method 
which  has  achieved  the  result  of  bit  rate  below  0.4,  and  mean- 
square  error  below  0.002. 

(2)  This  method  is  quite  simple  and  is  very  economical  for 
machine  cost  [1]. 

A  brief  introduction  to  recursive  IDPCM  is  shown  below: 

(1)  Take  coarse  subsample  spacing:  N  (e.g.,  N  =  8)  ; 
quantize  subsamples. 


(2)  Calculate  the  interpolated  value  of  intermediate  point 
between  two  subsamples;  use  this  interpolated  value  to  calculate 
the  difference,  and  quantize  the  difference. 

(3)  Use  the  intermediate  interpolated  value  added  to  the 
quantized  difference  to  calcualte  a  pair  of  intermediate  inter¬ 
polated  values  with  its  neighbor  two  subsamples,  etc. 

A  mode  is  defined  as:  S,  01,  02,  D3,  where  S  is  the  number 
of  bits  for  quantizing  subsamples,  and  DI,  D2,  and  D3  are  numbers 
of  bits  for  quantizing  differences.  But  01,  D2,  and  03  have  dif¬ 
ferent  numbers  of  bits.  For  example,  take  coarse  subsample 
spacing  N=8,  as  in  the  following: 


t 

bits 
-D^  bits- 


8x8  sub image 

0  —  subsamples  using  S  bits  quantization 
®  —  the  first  set  quantized  differences  using  bits 
+>  —  the  second  set  quantized  differences  using  D2  bits 
•  —  the  third  set  quantized  differences  using  bits 

Fi gure  1 . 


As  explained  above,  this  method  is  used  to  compress  an  entire 
image  without  regard  to  the  amount  of  detail  in  any  particular 
area.  Generally,  finer  sampling  is  required  in  the  neighborhoods 
of  sharp  gray-level  transitions,  while  for  relatively  smooth 
regions,  coarse  sampling  is  acceptable.  Although  recursive  IDPCM 
has  a  higher  data  compression  ratio  and  low  mean-square  error 
than  other  methods,  there  is  still  room  for  improvement. 

One  improvement  to  recursive  IDPCM  is  to  use  an  adaptive 
scheme.  With  this  method,  an  image  can  be  divided  into  subimages, 
where  a  high  bit  rate  is  needed  to  deal  with  relatively  complex 
subimages  but  a  lower  bit  rate  is  sufficient  for  relatively 
simple  subimages.  Here,  the  complex  subimages  are  defined  as 
neighborhoods  of  sharp  gray-level  transitions  and  simple  sub¬ 
images  as  neighborhoods  of  smooth  gray-level  transitions. 

Identification  of  Subimages 

In  order  to  use  an  adaptive  scheme,  we  have  to  detect  the 

complexity  level  of  each  subimage.  "One  way  to  measure  the 

o 

redundancy  of  an  image  and  to  compare  it  to  the  nominal  N  p  (NXN 
is  the  size  of  an  image,  p  is  bits  of  per  pixel)  bits  is  the  use 
of  the  histogram  statistics  and  the  associated  entropy  stati sties" [3] 
Entropy  represents  the  amount  of  information  associated  with  the 
set  of  coder  input  values  and  gives  a  lower  bound  on  the  average 
number  of  bits  required  to  code  those  inputs.  If  the  set  of  coder 
input  levels  is  Wl,  W2,  W3,...Wm  with  probabilities  pi,  p2,  p3,...pm, 
then  it  is  not  possible  to  code  them  without  distortion  using  less 


However,  this  requires  a  very  large  amount  of  computation. 
Some  approximation  has  to  be  introduced. 

Measuring  the  level  of  sharp  gray-level  transitions  can  be 
done  by  calculating  the  function: 


where  p^  is  the  mean  intensity  of  an  NxN  size  subimage,  and  pij 
is  the  intensity  of  each  pixel  in  a  subimage.  But  for  subimages 
a,  b,  and  c  shown  in  fig.  2,  we  get  the  same  value  for  R. 
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In  fact,  subimages  a  and  b  are  in  simple  categories  (low 
frequency)  and  the  subimage  c  is  in  complex  category  (high  fre¬ 
quency).  The  more  simple  an  image  is,  the  smaller  the  differences 
are  between  the  image  and  a  value  which  interpolates  it.  There¬ 
fore,  complexity  could  be  associated  with  the  differences.  The 
f uncti on , 


R'  = 


m 

A  2 

Z 

(  Pi 

-  p.  >2 

i 

M 

can  be  used  to  approxi mately  measure  the  complexity  of  a  subimage, 

A 

where,  p..  is  an  original  pixel  and  pi  is  an  interpolator  of  the 
subimage.  M  is  the  number  of  differences  in  a  subimage.  In 
recursive  IDPCM,  the  middle  interpolated  image  between  two  sub¬ 
samples  is  subtracted  from  the  original  pixel  and  the  first  set 
of  differences  are  quantized  at  bits.  Then  the  middle  inter¬ 
polator  added  to  its  quantized  difference  and  the  neighboring 
subsample  S  are  used  to  calculate  the  second  set  of  interpolated 
values  and  their  differences.  Therefore,  several  quantization 
levels  can  be  used  both  for  the  first  set  differences  and  second 
set  differences.  Fig.  10,  11,  and  12  show  that,  differences,  as 
defined  above,  can  represent  the  complexity  of  subimages.  This 
is  an  effective  and  feasible  method,  which  is  easy  to  compute. 


Quanti zati on 

"A  quantizer  is  a  device  whose  output  can  have  only  a  limited 


number  of  possible  values.  Each  input  is  forced  to  one  of  the 
allowable  output  values.  One  way  to  accomplish  this  is  to  divide 


Figure  5 


In  Recursive  IDPCM,  the  typical  distribution  of  differences 


is  shown  in  Figure  6, 


jP(x) 

I 


P(  x  )  is  the  probability  density  function  of  differences 
falling  into  A  X  range.  X  is  the  value  of  a  difference. 


the  input  range  into  a  number  of  bins  as  illustrated  in  Fig.  5. 

If  an  input  falls  into  the  Kth  bin,  the  output  is  the  value 
corresponding  to  the  center  of  Ith  bin  so  that  each  input  is 
rounded  off  to  the  center  of  the  bin  into  which  it  falls.  A 
uniform  quantizer  is  one  in  which  all  bin  widths  are  equal.  Non- 
uniform  quantizers  allow  different  bins  to  have  different  widths." 

Let  X  represent  any  input  value,  and  let  Q..  be  the  corresponding 
output  of  the  quantizer.  "If  all  values  of  X  within  the  bins  are 
not  equally  likely  then  the  squared  error  (X.  -  Q.)  must  be 

'  J 

weighted  by  the  probability  density  function  p ( X ^ )  "[5]. 

F  (X.  £  X  X.  +  d  X  ) 

P(  X.  )  *  Lim  - - 

1  4X->0  Z*  X 

F(Xi)  represents  the  number  of  differences  which  have  values 
between  [X^ ,  Xi  +  AX].  If  we  choose  AX  3  1  (because  of  digital 
signal),  pCX^)  is  approximately  equal  to  the  probability  density 
f uncti on . 

In  adaptive  IDPCM  the  quantization  strategy  is  to  choose  the 
quantized  levels  (Q^)  so  that  they  minimize  the  total  quantizer 
mean-square  error.  This  error  is  defined  as 

n  /-x.+l  - 

e*  *  2  <  x  -  Qi  >  P(X)  dx  (3) 

i-1  -'xi 

Taking  a  partial  derivate  of  equation  (3)  with  respect  to 


Q.  gives 
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Since  we  have  two  unknown  variables  Q^and  ,  we 
also  take  the  partial  derivative  of  equation  (3)  with  respect 
to  X± 


d  xi 


a  x.  JxL- 1 


(x"Qi-i)  p(x)dx 


xi+i 


(X  -  Q  }  P  (X)  dX 


and  get  (Xj_  -  Q^) 2  P(X)  -  (X±  -  Q^2  P(X)  *  0 


Xi  "  Qi-1  3  ^  ’  Xi 


From  simultaneous  equations  (5)  and  (6) 
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optimal  quantizer  output  levels  are  achieved,  which  is  shown 
in  Fig-.  7 


Quantization  Input  Range- 


In  Fig.  7,  the  original  image  is  sampled  at  512x512  resolu¬ 
tion  with  8-bits  per  pixel.  The  image  data  compression  steps  were 
as  follows: 

(1)  The  original  image  was  subsampled  at  every  8th  line  and 
every  8th  pixel.  Each  subsample  was  quantized  to  6-bits  with  a 
uniform  quantization.  The  maximum  and  minimum  quantization  levels 
were  256  and  0. 

(2)  The  middle  point  interpolated  pixels  were  subtracted  from 
the  original  pixels,  and  the  differences  were  optimally  quantized 
in  3  bits  if  the  mean-square  differences  in  a  subimage  were  larger 
than  4.0,  (the  average  error  in  subsamples  caused  by  quantizing 
8-bit  original  pixels  to  6-bit  subsamples).  Of  course,  those 
differer.  es  were  not  stored  if  the  mean-square  differences  in  a 
subimage  were  smaller  than  4.0. 

(3)  The  subsamples  and  the  quantized  middle  point  differences 
were  then  used  to  calculate  the  second  set  differences.  Suitable 
thresholds  for  the  mean-square  differences  in  a  subimage  were 
chosen  to  determine  how  many  bits  are  being  used  for  quantizations 
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detail  subimages. 


(4)  The  subsamples  and  quantized  differences  according  to  the 
R'  were  used  to  reconstruch  the  image.  Fig.  8  shows  the  original 


image.  Fig.  9  shows  the  processed  image  by  using  recursive  IDPCM, 
and  Fig.  10  shows  the  image  by  using  adaptive  recursive  IDPCM. 
Comparing  Fig.  9  with  Fig.  10,  we  can  see  that  the  details  in 
Fig.  10  did  improve. 

The  total  bit  requirements  for  the  images  is  the  sum  of  bits 
for  subsamples,  differences,  and  the  modes  of  each  subimage. 

Since  the  number  of  bits  for  representing  modes  of  each  subimage 
is  very  small,  these  extra  bits  can  be  negligible, 
total  bits  =  3x3cN-j  +  6XN2  +  2x2xN^  +  3xl2xN^ 

N 1  =  number  of  subimages  of  R‘  >  4.0,^  =  number  of  total  sub¬ 
images,  N ^  =  number  of  subimages  of  its  R'  (the  second  set  dif¬ 
ferences)  between  4.0  and  12.0,  3  number  of  subimages  of  R' 

(the  second  set  differences)  >12.0.  In  adaptive  recursive  IDPCM 
for  the  image  shown: 

Bit  Rate  of  Per  Pixel  =  0.35259 
Mean-Square  Error  =  0.001798 
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( V. )  Applications  of  Stationary  Transform  Processin 


Space-Variant  Image  Processing 

This  work  has  focused  on  three  areas:  (1)  Improved  estimates 
of  local  autocovariance  statistics  for  direct  implementation  of 
spati al ly-adapti ve  image  processing  algorithms.  (2)  Development 
of  adaptive  image  restoration  algorithms  for  nonstationary  images 
degraded  by  blur  and  additive  noise.  (3)  Development  of  adaptive 
image  enhancement  algorithms  for  local  contrast  transformation, 
with  application  to  color  imagery. 

The  results  of  this  work  are  reported  in  the  sections  to 
follow.  In  addition.  Appendix  A  gives  details  of  an  extension 
to  the  geometri c- transform  algorithm  (see  1982  and  1983  AFOSR 
Annual  Reports).  This  computer  algorithm  is  capable  of  generating 
a  grid  whose  cells  are  of  specified  area  and  unspecified  shape. 

The  procedure  is  iterative,  producing  a  better  solution  after  each 
iteration,  and  converging  even  for  quite  large  differences  in 
cell  area. 

Improved  Estimates  of  Local  Autocovariance  Functions 

Our  previous  work  in  this  area  used  the  first-order  Markov 
model  for  autocovariance  functions  estimated  from  image  subblocks. 
First-order  parameters  were  obtained  by  fitting  an  exponential 
curve  to  the  autocovariance  data.  Recent  tests,  however,  indicate 
that  this  approach  gives  highly  variant  statistical  estimates,  due 
to  the  inclusion  of  unreliable  high-order  lag  data  in  the  curve 
fitting  process. 


In  Appendix  B  we  report  results  from  our  experiments  to 
find  useful  measures  of  local  image  autocovariance  parameters 
from  small  subblocks  of  data.  The  basis  of  this  work  is  a  new 
technique  for  estimating  the  correlation  parameters  of  first- 
order  Markov  (nonseparable)  exponential  autocovariance  models. 

The  method  assumes  that  image  data  is  stationary  within  NxN-pixel 
subblocks.  A  measure  of  the  usefulness  of  the  ensuing  correlation 
parameters  may  be  had  by  observing  their  correlation  with  signal 
activity  within  the  scene.  Subblocks  of  dimension  N*16  are 
shown  to  provide  estimates  fitting  this  criteria,  even  when  the 
original  data  is  degraded  by  significant  amounts  of  blur  and 
noise.  This  latter  fact  bodes  well  for  spati al ly-adapti ve  image 
restoration  applications.  An  application  to  bl ock- i nterpol ati ve 
data  compression  is  also  included  in  Appendix  B. 

Adaptive  Image  Restoration 

The  basic  concept  of  stationary  transforms  in  image  restora¬ 
tion  is  as  follows: 

.  Assuming  the  image  to  be  degraded  by  blur  and  additive 
noise,  then  - 

.  Estimate  the  autocovariance  parameters  by  the  method 
outlined  previously. 

.  Apply  geometric  transforms  to  produce  stationarity  in 
correlation  length. 

.  Filter  the  resulting  data  with,  e.g.,  a  Wiener  filter 
based  on  the  stationary  autocovariance  distribution. 


-.y.y.y 


,\*V 


.  Apply  the  inverse  geometric  transform  to  restore  the  original 
geometry. 

This  approach  is  quite  feasible  for  1-D  data  processing,  sine 
the  geometric  transforms  simply  translate  into  straightforward 
decimation  and  interpolation  operations.  However,  we  showed  in 
the  1983  AFOSR  Annual  Report  that  the  2-D  implementation  of  the 
above  steps  results  in  highly  distracting  artifacts  in  the  form 
of  regular  patterns  of  correlated  noise.  In  an  attempt  to  over¬ 
come  this  limitation  we  investigated  a  modification  of  the  above 
technique : 

.  Estimate  autocovariance  parameters  in  NxN-pixel  subblocks 
as  before. 

.  Subsample  (decimate)  subblocks  by  a  linear  scaling  factor 
proportional  to  these  parameters  in  order  to  produce 
constancy  over  all  the  parameter  set. 

.  Place  the  new  subblocks  at  the  centers  of  the  original 

subblocks,  and  apply  a  nonadaptive  Wiener  filter  as  before. 

.  Restore  the  subblocks  to  their  original  dimensions  by 
i nterpol ati on . 

Figures  1(a)  and  1(b)  show  the  results  of  this  procedure  for 
the  case  where  additive  noise  is  the  only  source  of  degradation. 
The  resulting  restoration  in  Figure  1(b)  is  encouraging,  since 
the  residual  blocking  effect  can  most  likely  be  removed  by  band¬ 
pass  filtering  incorporated  at  the  same  time  as  the  noi se- f i 1 ter- 
ing. 

We  have  also  investigated  several  new  approaches  for  adaptive 
image  restoration  using  space- varyi ng  filters.  A  simple  type  of 


space-vari ant  deblurring  is  implemented  via  an  Adaptive  Biased 
Laplacian  operator.  This  is  a  3x3  kernel  given  by: 


L(x,y) 


-1 

-1 

-1 

-1 

8+d(x,y) 

-1 

-1 

-1 

-1 

(1) 


The  deblurring  is  given  by  a  convolution  between  L(x,y)  and 
the  image  i(x,y).  The  term  d(x,y)  is  a  space-varyi ng  bias  which 
adapts  to  local  image  statistics;  for  example,  image  variance 
var(x,y)- 

d( x  ,y )  =  f ( var(x,y) )  (2) 

where  f(  )  is  moni toni cal ly  decreasing. 

When  var(x,y)  is  high,  for  example  at  large  edges  we  want 
L(x,y)  to  be  a  high-pass  filter.  In  regions  of  low  variance,  we 
want  L9x,y)  to  be  all-pass.  These  constraints  are  satisfied  when 

d(x,y)  =  1  +  9xexp(-fxvar(x,y)2)  ;  f>0  (3) 

From  eqn.  (3),  we  notice  that  the  minimum  value  of  d(x,y)  is 
unity,  and  the  maximum  is  10.  In  conjunction  with  f=0.011,  this 
range  is  found  to  give  good  results. 

Normalization  of  L(x,y)  is  necessary  to  avoid  large  spatial 
gain  variations  during  the  convolution:  the  resulting  normalized 
kernel  is  given  by; 

L 1 ( x ,y )  =  L(x,y)  v  d(x,y)  (4) 

The  adaptively-processed  image  is, 

v  •  ■  ^ 


g(x,y)  =  i ( x ,y )  **  L' (x,y) 


(5) 


Separating  terms  in  equ.  (5),  we  have; 
g(x,y)  =  i (x,y)  **  (L/d(x,y)  +  i(x,y)  (6) 

where  L  is  the  usual  3x3  Laplacian  kernel  consisting  of  -1  in  all 
elements  except  the  center  element  of  8.  An  additional  weight 
in  the  first  term  prevents  saturation,  i.e., 

9(x,y)  =  0 . 3[i ( x ,y )  **  (L/d(x,y))]  +  i(x,y)  (7) 

Figures  2(a)  and  2(b)  show  an  image  before  and  after  pro¬ 
cessing  with  equ.  (7).  The  effect  of  unwanted  noise  boosting  in 
areas  of  low  signal  activity  is  bypassed  in  this  adaptive 
approach . 


Adaptive  Imaqe  Enhancement 


A  new  algorithm  for  adaptive  contrast  enhancement  is  described 
in  this  section.  The  aim  of  the  algorithm  is  to  compensate  for 
nonuniform  illumination  in  scenes,  and  is  based  on  the  often 
observed  fact  that  local  maps  of  mean  and  variance  correlate 
quite  strongly  in  the  positive  sense.  The  philosophy  of  the 
algorithm  may  be  outlined  as  follows: 

.  If  the  local  mean  T(x,y)  is  high,  then  leave  the  image 
unal  tered. 

.  If  the  local  mean  T(x,y)  is  low,  then  increase  the  variations 
about  the  mean  as  well  as  the  mean  itself. 

For  convenience  we  will  drop  the  (x,y)  image  coordinates  and 
represent  the  algorithm  as  follows: 

g  =  k ^ ( i  -  T)  +  k2  +  T  (8) 

The  terms  k,  and  k,  are  functions  of  T(x,y).  In  order  to  satisfy 
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the  algorithm  specifications  they  should  exist  within  the  following 


1  imi  ts 


k1  =  <  1  ;  i  =  255 

<  3  ;  T  =  0 
max(k^)  s  3 

k2  “  <  0  ;  T  =  255 

<  C  ;  T  =  0 


Functional  forms  which  accomodate  these  limits  are  given  in  a 


general  form  by. 


kl  =  AC  "2^5)  +  1 

1  t 

k2  =  C(  1  -255 ; 


Optimum  values  for  the  constants  based  on  experimental  trials 
are:  A  =  100,  C  =  70,  t  =  2 . 

The  terms  k^  and  k2  are,  respectively,  the  local  contrast 
and  local  mean  adjustment  factors.  Calculation  of  local  means 
T  is  achieved  by  neighborhood  averaging  over  9x9-pixel  windows, 
Figures  3(a)  &  (b)  illustrate  the  effect  of  this  adaptive 
image  enhancement.  Individual  red,  green,  and  blue  bands  were 
processed  separately  to  produce  the  color  imagery. 


Note 


Appendix  A  was  submitted  for  publication  in  "Computers  and 
Graphics",  January  1984. 

Appendix  B  was  submitted  for  publication  in  IEEE  Trans,  on 
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I.  INTRODUCTION. 


This  paper  addresses  a  simply-posed  problem:  starting  with  a  simple 
grid  of  square  cells,  determine  the  coordinates  of  a  new  grid  in  which  the 
area  of  every  cell  is  specified  individually.  The  shapes  of  cells  in  the  new 
grid  are  unspecified. 

This  problem  arose  in  an  image  processing/ remote  sensing  application  [1]. 
We  envisage  that  the  solution  may  find  other  applications;  for  example,  in  com¬ 
puter  graphics,  mapping,  integrated  circuit  fabrication... 

The  flavor  of  the  paper  is  heuristic.  Convergence  of  the  algorithm  is 
demonstrated  by  examples,  without  proof. 

II.  ALGORITHM. 

Figure  1  illustrates  the  broad  objectives  and  constraints  of  the  problem. 
Each  unit-area  cell  in  the  original  regular  grid  is  to  be  expanded  to  form  an 
arbitrary  shaped  quadrilateral  of  specified  area.  No  other  constraints  are 
imposed,  apart  from  the  obvious  ones  of  contiguity  and  connectivity  in  the 
transformed  grid.  Simply  stated,  we  wish  to  know  the  coordinates  of  each 
cell  in  the  transformed  grid,  relative  to  an  arbitrary  set  of  orthogonal  axes 
in  the  x-y  plane.  If  we  define  a  regular  grid  of  control  points  whose  coor¬ 
dinates  lie  at  the  vertices  of  each  cell,  then  the  problem  is  to  determine  how 
such  a  regular  grid  is  distorted  by  the  many  interactions  of  local  area  ex¬ 
pansions  centered  on  each  cell. 

Our  approach  is  to  superpose  the  expansion  of  each  cell  on  the  global 
grid,  subject  to  the  constraint  that  each  expansion  should  minimally  affect 
the  relative  geometry  of  other  cells.  This  philosophy  is  easily  explained  by 


a  1-D  illustration.  In  Figure  2(a)  we  begin  with  uniformly  spaced  control 
points.  In  Figure  2(b)  an  expansion  factor  of  2  is  applied  to  the  center 
of  a  cell  at  xQ,  resulting  in  localized  stretching.  The  next  expansion 
is  applied  to  the  neighboring  cell,  and  so  on.  Continuing  this  process 
yields  a  1-D  control  point  array  whose  coordinates  reflect  the  desired  ex¬ 
pansion  of  cells.  No  iterations  are  required.  The  basic  1-D  geometric 
transformation  can  be  represented  by  a  spatial  displacement  ax. 

Ax  s  (E  -  1)  (la) 


where  E  is  an  "expansion  factor,"  and 

Compression  -  0  s  E  <  1 
Expansion  -  1  s  E  <  • 


(lb) 


In  the  2-D  case,  it  is  not  possible  to  design  a  transformation  which 
expands  the  area  of  a  single  cell,  while  preserving  the  geometry  of  all 
other  cells.  What  is  needed  is  a  technique  for  reducing  the  mutual  interference 
of  local  cell  expansions.  Consider  a  single  radial  expansion  centered  on  a 
single  cell  in  a  grid  of  control  points,  as  shown  in  Figure  3.  Straight¬ 
forward  radial  expansion  takes  place  within  a  circle  of  radius  r  centered 
on  the  expansion  cell.  Elsewhere,  more  remote  control  points  are  simply  dis¬ 
placed.  Hence,  we  can  define  a  geometric  transformation  as  follows: 


AX 

Ay 


x(E-l) 

y(E-l) 


(2a) 


-2- 


where  x,  y,  and  e  are  control  point  coordinates  relative  to  the  center  of  ex- 
pansion  (see  Figure  3).  Nominally,  r  =  d//?  in  order  to  confine  expansion 
to  a  single  cell  as  much  as  possible.  The  result  of  applying  a  linear  expan¬ 
sion  of  E  a  2  to  a  single  cell  is  shown  in  Figure  4.  Note  that  the  area 

2 

expansion  of  the  cell  is  E  ,  or  4.  Also,  notice  that  the  influence  of  this 
expansion  on  other  cells  diminishes  at  greater  distances  from  the  expansion 
center. 

We  have  found  that  the  simple  geometric  transformation  in  Eqn.  (2)  works 
well  if  applied  sequentially  to  all  cells  in  the  grid,  provided  that  correction 
factors  are  used  to  compensate  for  the  influence  of  single-cell  expansions 
on  all  other  cells.  Furthermore,  iterating  the  procedure  yields  the  desired 
cell  ar»as  in  most  cases. 

The  iterative  algorithm  proceeds  as  follows: 

Notation 

Assume  initially  an  n  x  m  grid  of  unit-area  cells. 

C^.  -  cell  (i,j)  (1  s  i  s  n  :  1  s  j  s  m) 

E^j  -  desired  linear  expansion  of  C.j 

Ejj  -  desired  area  of 

A..  -  actual  area  of  C. .  at  any  time. 

'  J  ■  J 


Algorithm: 


rSTART 


For  1  =  IORDER  (1)  IORDER  (n) 

For  j  =  JORDER  (1)  -*•  JORDER  (m) 
Measure  A^. 

Form  the  corrected  expansion  E  =  E. 
Apply  Eqn.  (2)  to  C..,  using  E. 

*  w 


L-  END  OF  LOOPS 


RETURN  TO  START,  repeat  for  K  iterations. 


The  arrays  IORDER  and  JORDER  are  needed  to  generalize  the  order  in  which 
the  cells  are  treated,  in  other  words,  each  cell  is  considered  only  once 
during  an  iteration,  but  cells  may  be  addressed  in  any  desired  sequence.  We 
will  refer  to  this  as  the  cell  sequence. 

NOTE:  In  applying  Eqn.  (2),  we  measure  x  and  y  from  the  specified  control 
point  to  the  centroid  of  C^. 


III.  EXAMPLES* 

Most  of  our  examples  concern  a  grid  of  16  cells;  i.e. ,  25  control  points. 
The  arbitrary  set  of  expansions  shown  in  Figure  6  is  used  throughout  the  tests. 
Several  variations  of  the  iterative  algorithm  are  presented  in  order  to  demon¬ 
strate  the  following: 

•  effect  of  cell  sequence. 

•  effect  of  radius  r  in  Eqn.  (2). 

•  existence  of  an  upper  limit  on  E^j. 

*  The  figures  in  this  section  are  not  drawn  to  scale. 


In  all  tests,  we  employ  an  rms  expansion  error  measure  for  the  solution, 
given  by: 

(l  n  m  ?  0)l/2 

±  A  sh  tEu  -  V/  <3> 

where  k  is  the  iteration  index. 

Cell  Sequence. 

Figure  5  shows  the  notation  for  cell  sequence.  Tests  indicate  that  row- 
by-row  or  column-by-column  sequences  are  superior  to,  say,  circular  or  spiral 
ordering.  Indeed,  for  the  particular  16-cell  example  selected,  only  the 
former  yield  convergence  in  the  solutions.  Figures  6(a)  and  (b)  show 
the  transformed  cells  after  100  iterations  using  the  following  cell  sequences, 
(i)  Fig.  6(a)  -  1  2  3  4,  5  6  7  8,  9  10  11  12,  13  14  15  16. 

(ii)  Fig.  6(b)  -  1  2  3  4,  8  7  6  5,  9  10  11  12,  16  15  14  13. 

In  both  (i)  and  (ii),  the  nominal  value  r  *  6//Z~  was  used.  Evidently,  the 

rate  of  convergence  and  form  of  solution  is  affected  by  the  cell  sequence, 
although  only  slightly  in  this  example.  However,  these  effects  depend  greatly 
on  the  expansion  values.  In  some  cases,  even  the  sequence  in  (ii)  does  not 
yield  a  solution. 

Radius  parameter,  r. 

The  nominal  value  of  r  *  d//Z~  is  not  optimum.  For  most  expansion  sets, 
a  value  can  be  found  which  dramatically  increases  the  convergence  rate.  For 
example,  a  value  of  r  =  /2d  in  conjunction  with  the  cell  sequence  in  (ii) 
gives  a  solution  with  e ( 1 00)  *  7  x  10"®;  i.e. ,  a  factor  of  100  better  than 
the  nominal  r  solution.  The  resulting  grid  of  cells  is  shown  in  Figure  6(c). 
Clearly,  this  value  of  r  also  leads  to  a  somewhat  different  solution. 


Upper  limit  on  expansions. 

The  series  of  results  in  Figures  6(d), (*),  and  (f)  illustrates  the 
behavior  of  the  algorithm  as  the  expansion  applied  to  one  particular  cell 
is  increased.  Ultimately,  the  boundaries  of  adjacent  cells  cross,  creating 
additional  cells  (cell  nos.  15  and  16  in  the  example  shown).  Our  experience 
is  that  the  maximum  allowable  expansion(s)  is  a  function  of  r  anc[  the  cell 
sequence. 

Some  additional  interesting  results  are  now  discussed.  Tests  of  the 
iterative  algorithm  using  symmetrical  expansion  sets  yield  Figures  '7(a) 
and  (b).  We  observe  that  the  asymmetry  of  the  cell  sequence  is  responsible 
for  the  asymmetry  of  the  solutions.  An  example  using  64  cells  is  given  in 
Figure  0. 


IV.  REMARKS. 

We  have  demonstrated  the  behavior  of  an  iterative  procedure  for  generating 
connected  quadrilaterals  of  specified  area.  The  solutions  are  nonunique,  and 
depend  on  two  free  parameters  -  radius  r  in  Eqn.  (2),  and  the  cell  sequence. 
Obviously,  the  nonuniqueness  is  a  consequence  of  the  overdeterminacy  of  the 
problem.  Additional  constraints  on  the  quadrilateral  cells,  such  as 
shape  constraints,  would  reduce  the  number  of  solutions.  The  algorithm  pre¬ 
sented  here  usually  converges,  even  when  the  expansion  set  spans  a  large  dynamic 
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Result  using  8x8  grid  of  expansions 
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(  e(50)  =  0.005  ) 


I .  INTRODUCTION 


the  design  of  digital  image  processing  algorithms  is  critically  dependent 
an  the  spatial  statistics  of  the  imagery.  The  largest  class  of  algorithms 
assumes  stationary  statistics!  for  example,  the  Wiener  filter  for  least-squares 
image  restoration  Cl3,  and  the  DPCM  coder  in  data  compression  [2]  .  The 
main  advantage  of  such  algorithms  lies  in  the  relative  simplicity  of  their 
design  and  implementation.  A  drawback,  however,  is  that  images  axe  inherently 
nanstationary;  consequently,  the  results  of  processing  are  suboptimal ,  Hie 
Wiener  filter,  for  example,  is  designed  from  estimates  of  the  object  (and 
noise)  autocovaxlance  functions.  The  question  is:  how  do  we  obtain  meaningful 
estimates  of  the  abject  covariance  function?  Typical  imagery  consists  of 
large  regions  with  nearly  uniform  intensity,  plus  a  relatively  few  «w»n  areas 
of  high  signal  activity,  such  as  edges,  natural  texture,  etc.  Therefore, 
an  estimate  of  the  autocovariance  based  on  the  global  image  data  would  yield 
a  rather  broad  function,  since  most  of  the  data  tends  to  be  highly  correlated. 

A  Wiener  filter  designed  from  such  an  estimate  would  yield  a  very  smooth 
restoration,  and  important  edge  detail  would  remain  blurred.  Since  the 
Wiener  filter  provides  a  compromise  between  noise  smoothing  and  deblurring, 
with  a  fairly  conservative  bias  towards  noise  smoothing,  a  better  approach 
is  to  estimate  autocovariance  from  regions  of  low  correlation,  i.e.  edges, 
etc.  The  resulting  restored  images  exhibit  a  greater  degree  of  sharpness, 
althou^x  they  are  inevitably  noisier  than  normal  least-squares  restorations. 
Nevertheless,  *11  of  these  techniques  are  handicapped  because  they  do  not 
account  for  nanstationarity  in  the  image  data. 

A  more  recent  class  of  image  processing  algorithms  accommodates  this 
nanstationarity.  These  algorithms  share  the  characteristic  of  adapting  to, 
or  tracking,  the  local  variations  in  image  correlation.  The  resulting  pro¬ 
cessors  are  space-varying ,  in  accordance  with  the  image  nanstationarity. 


Lahart  [3]  has  reported  a  method  of  implementing  local  least-squares  resto¬ 
ration  in  which  two  different  autocovariance  functions  were  used,  depending 
on  whether  pixels  were  considered  to  belong  to  high  or  low  signal -activity 
categories.  Anderson  and  Ne travail  [4]  used  a  masking  function  to  implement 
a  form  of  space-varying  noise-smoothing  which  adapted  to  local  signal  details. 
A  different  approach  to  achieving  the  same  goal  is  due  to  Kasturi  CS],  who 
considered  the  particular  case  of  signal -dependent  noise.  (See  also  Fzoehlich, 
Walk  up ,  and  Asher  [6].)  The  benefits  of  space-varying  processing  have  proved 
even  more  dramatic  in  the  field  of  image  data  compression.  A  review  of  adapt¬ 
ive  coding  methods  is  provided  by  Habibi  [7],  All  of  the  work  cited  above 
depends  on  some  kind  of  direct  estimation  of  local  image  statistics.  Other 
sethods  use  nonstationary  models  indirectly}  an  excellent  example  is  Widrow's 
IMS  adaptive  filter  C  3] . 

In  order  to  quantify  the  nonstationarity  of  images,  and,  often,  in  order 
to  implenent  adaptive  processing,  we  require  estimates  or  measurements  of 
ths  local  image  statistics,  specifically  the  auto  covariance  function.  The 
simplest  way  of  achieving  this  is  to  divide  the  image  into  N  x  N  -  pixel 
subblocks  (»  “  16)  ,  an d  calculate  the  standard  biased  or  unbiased  autocovar¬ 
iance  function  of  each  subblock.  In  affect,  each  subblock  is  treated  as  part 
of  awide-sense  stationary  field.  It  is  well-known,  however,  that  reliable 
power  spectral  estimation  requires  much  larger  amounts  of  data.  Neverthe¬ 
less,  as  our  work  shows,  it  is  possible  to  obtain  useful  maps  of  local  auto¬ 
covariance  parameters  if  we  assisa e  simple  parametric  autocovariance  models. 
Specifically,  we  esplgy  the  popular  first-order  Markov  random  field  model 
for  data.  Three  topics  are  addressed  herein: 

•  estimation  of  local  autocovariance  parameters. 


optimum  choice  of  N. 


correlation  between  the  estimated  autocovariance 
parameters  and  observed  signal  activity. 

In  the  sections  to  follow,  we  review  the  autocovariance  model  and  discuss 
possible  ways  of  estimating  the  parameters  of  each  local  subblock.  Examples 
are  given  which  illustrate  seme  of  the  problems  encountered  when  the  estimates 
axe  applied  to  real  image  data.  Finally,  we  give  an  example  of  the  value  of 
local  estimates  in  adaptive  data  coapression. 


XX.  LOCAL  IMAGE  MODELS. 

Global  wi de-sense  stationary  image  data  is  generally  modeled  by  a  first- 
order  Markov  random  field  with  a  nans  ep  ax  able  exponential  autocovariance  func¬ 
tion  C9]  given  by 

C(k,£)  -  a2  exp  {  -  /~p  k*+  p  Z *  }  (1) 

1  2 


where  o2 
are  p^  » 


is  the  variance.  Typical  values  adopted  for  the  correlation  parameters 
p  ■*  0.0025,  resulting  in  the  simpler  isotropic  form 


C(k,f) 


a2  0.95 


/k2  +  Z2 


U) 


This  model  is  often  assigned  in  nen  adaptive  DP  CM  coding  systems. 

A  new  model  for  nenstationazy  image  data  f(i,j)  is  given  by: 

C(k,f:i,j)  -  o2(i,j)  exp  {  -  /p  U,j)k*  +  p  (i,j)  f4  }  (3) 


It  is  assumed  that  f(i,j)  is  stationary  over  any  N  x  N  -  pixel  subblock  cen¬ 


tered  on  the  coordinates  (i,j)  ,  and  that  values  of  o2,  p  ,  and  p  may  be 
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assigned  to  every  location  (i,j)  .  (In  practice ,  these  values  will  be 

correlated  with  their  neighbors .  Later  in  the  paper  we  will  discuss  the 

implications  of  this.)  The  next  section  discusses  procedures  for  estimating 

the  spatially-variant  autocovariance  parameters  p  (i,  j)  and  p  (i,j)  . 
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III.  ESTIMATION  OF  LOCAL  CORRELATION  PARAMETERS. 

The  starting  point  for  our  estimates  of  p  values  is  the  usual  biased 
autocorrelation  estimate, 

i  +  Hh  -  w  -i  j  +-f--  \t\  -i  (4) 

f(m,n)  f  (m+k,  n+£) 


R(k^si,j)  - 


a  ■  i  - 


N 


n  -  j  - 


N 


or  the  biased  autocovariance  estimate,  given  by, 

i  +  j+Hr  -  kl-  i 


C(k,Z:i,j)  « 


N 


a  ■  i  - 


N 


n  -  j  - 


N 


[f(m,n)  -  f  (i,j)l 
(5) 


Cf  (a  +  k,  n  +  l)  -  f  (i,  j)  3 


where  f  (i,j)  is  the  spatial  mean  of  an  N  x  N  -  pixel  subblock  centered  on  (i, j) 


Simple  parametric  estimates 


From  Eqn.  (3),  simple  parametric  estimates  of  and  p^  are  given  by 


p  (i.j)  •  { 


TTT  ^ 


r  ct^.osi,  jH « 


;  hen  k^  •  0  and  k^  •  1,  these  estimates  are  equivalent  to  the  first-order 
autoregressive  (AR)  estimates  found  in  Box- Jenkins  time  series  analysis  [10  J. 
For  example,  p  is  then 


P  (i» j)  m  l  In 
1 


p  C<0,0:i, j}~| 

L  C(l,0:i,j)J 


In  practice,  our  constraint  that  N  should  be  small  results  in  autocovariance 
data  with  a  high  degree  of  variance,  causing  the  p-estimates  to  be  inaccurate. 

A  A 

Furthermore,  the  values  of  C(l,0:i,j)  and  C(0,l:i,j)  can  be  (and  often  are) 
negative,  rendering  Eqns.  (6)  and  (7)  indeterminate.  For  example.  Figure  2 
shows-  a  map  of  p^(i,j)  obtained  by  applying  Eqn.  (6)  to  the  image  in  Fig¬ 
ure  1  with  k.  ■  1  and  k_  *  3.  The  map  was  obtained  by  computing  p  at  every 
*  2  1 

pixel  coordinate  in  the  original  image,  with  overlapping  subblocks  of  16  x  16 

pixels.  The  large  dark  areas  correspond  to  negativity  in  the  argument  of  £n 

in  Eqn.  (6),  caused  by  negative  values  of  either  C(k_,0:i,j)  or  Cfk_,0ti,j). 

±  2 

At  this  point  we  should  note  the  influence  of  the  choice  of  lags  k^  and 
k2  on  the  estimation  quality.  Estimation  theory  tells  us  that  the  lowest 
order  lags  provide  the  most  reliable  autocovariance  data  (least  variance)  [11  ] 


suggesting  that  k^  *  0,  k^  ■  1  be  the  optimum  data  to  use.  However,  the 
autovariance  at  zero  lag  o£ten  contains  a  significant  "error"  equal  to  the 
variance  of  unwanted  film-grain  noise  (assuming  the  image  data  originates 
from  film).  This  explains  the  choice  of  k^  above.  The  selection  of  k^  in¬ 
volves  a  tradeoff  between  the  greater  reliability  of  low  order  lag  data  and 
the  increased  error  sensitivity  of  Sqn.  (€)  when  is  small.  However,  this 
form  of  estimation  gives  poor  results  for  any  combination  of  lag  data,  as 
Figure  2  demonstrates. 


Modified  parametric  estimates. 


A  new  method  of  p -estimation  avoids  the  indeterminacy  problem  encountered 
in  the  simple  estimation  of  Eqn.  (6) .  For  discussion  purposes,  estimation  is 
executed  in  two  stages,  although  the  calculations  can  be  combined  in  a  single  step 

First,  the  image  is  radiometrically  transformed  [12],  [13]  to  yield  data 
exhibiting  block-stationary  mean  and  variance.  The  transform  is  given  by. 


9(i'j)  “  0(137  +  fa 


where  o(i,j,)  and  7(1, j)  are  measured  in  an  M  x  M  -  pixel  subblock  centered 
on  (i,j).  The  mean  and  variance  measured  over  any  M  x  M  -  pixel  suhblock 
of  g(i.j)  are  the  constants  7*  and  o|,  respectively.  Typically,  M 

is  8  pixels,  in  the  second  stage,  the  biased  autocorrelation  estimates 
in  Eqn.  (4)  are  used  to  form  the  ratios. 


zl(Jcl'lt2:i'^) 


*2(W1'j) 


ftCk^Osi.j) 

&{k2;osi,j) 


If  we  accept  the  exponential  autocovariance  model  in  Eqn.  (1),  then 


the  expected  value  of  z^  (say)  is  [14]  , 


2  *1 

Ca*  exp  >  ♦  (fg)  3  Cl  -  ] 

[oj  exp  {-^/"p^  J  *  (?.)2jCl  ~  -jj2*  ] 


(10) 


The  expected  value  of  z  is  likewise  related  to  p  .  Eqn.  (10)  represents 
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a  nonlinear  relationship  between  z^  and  /  p  which  can  be  implemented  in 
the  fora  of  a  look-up  table  for  the  purpose  of  obtaining  p  values.  (See 
Figure  3.) 

Figures  4(a)  and  4(b)  show  the  distributions  of  z^  and  z^  obtained  by 

applying  Eqn.  (9)  to  the  autocorrelations  of  the  transformed  data  g(i,j). 

The  parameters  of  Eqs.  (4),  (8)  ,(9)  and  (10)  are:  N-16,  M  -  8,  a2  m  1500, 

s 

-  128,  m  1,  and  k2  -  3.  The  histogram  distributions  of  z^  and  z2 
values  are  depicted  in  Figure  3  relative  to  the  z  look-up  table.  Most 
(80%)  of  the  values  fall  between  the  limits  1.154  and  1.194  corresponding 
to  the  range  0.0  sfp  s  0.25.  This  is  probably  the  useful  range  of  p  for  ap¬ 
plications  in  image  processing  and  analysis.  Values  of  z  outside  this  range 
are  caused  by  breakdown- of  the  exponential  model  in  Eqn.  (3).  For  example, 
in  certain  subbloctathe  autocovariance  est-imate  in  Eqn.  (5)  contains  3trong 
periodic  components,  interestingly,  the  mode  of  the {p -histogram  occurs  at 

{p  «  0.1,  corresponding  to  the  correlation'  coefficient  exp  (-0.1)  «  0.9  often 
used  in  nonadaptive  processing. 

A  significant  advantage  of  this  new  p -estimator  is  that  z^  and  z2  are 
always  positive,  thereby  avoiding  the  problem  caused  by  negative  autocovari¬ 
ance  data  in  Eqn..  (6).  Furthermore,  the  estimates  are  much  less  sensitive 
to  very  small  values  of  autovariance  at  lag  *2 ,  caused  by  random  errors  in 
the  estimation  of  autocovariance  itself  (Eqn.  (5)).  Also,  estimates  may 
be  improved  by  pre-processing  the  maps  of  z  and  z,  using  a  smoothing  algor- 


ithra,  thereby  reducing  the  random  errors  incurred  in  the  model-fitting  process 

The  effect  of  neighborhood-averaging  is  evident  in  the  narrowing  of  the 

histograms  of  z^  and  z 2  also  shown  in  Figure  3.  (Values  of  z  outside  the 

useable  range  in  Figure  3  may  be  compressed  by  means  of  a  slightly  modified 

look-up  table.  This  form  of  biasing  to  achieve  reduced  variance  in  parameter 

estimates  is  common  in  spectral  estimation.)  Application  of  the  previous 

look-up  tables  (Figure  3)  produces  the  maps  of  /  p  and  /  p  in  Figures 
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5(a)  and  5(b).  An  important  feature  of  these  maps,  also  apparent  in  examples 
later  in  the  paper,  is  the  correlation  between  the  autocovariance  parameters 
and  signal  activity.  It  is  interesting  to  note  that  maps  of  local  variance  - 
often  employed  for  signal  measures  -  are  not  as  sensitive  to  low  contrast 
detail.  The  indication  is  that  the  nonseparable  exponential  model  for  auto¬ 
covariance  is  worthwhile  for  describing  nonstationary  data  if  applied  locally. 

IV.  BEHAVIOR  OF  THE  ESTIMATOR. 

This  section  examines  important  practical  aspects  of  the  behavior  of 
Eqn.  (9) . 

Subblock  dimension  (H) . 

Local  autocovariance  estimation  is  based  on  small  subblocks  of  N  x  N 
pixels,  each  one  of  which  is  assumed  to  belong  to  a  global  wide-sense  station¬ 
ary  random  field.  Attaining  useful  resolution  of  the  image  nonstation arity 
requires  that  N  be  small,  and  yet  this  violates  the  normal  rules  of  power 
spectral  estimation.  This  would  suggest,  then,  that  any  results  of  local 
estimators  are  useless.  However,  the  fact  that  our  estimates  correlate  well 
with  observed  signal  activity  suggests  the  contrary.  Further  support  for 
Eqp.  (9)  is  now  given. 


Consider  a  one-dimensional  window  of  data  which  contains  a  single  edge 

(see  Figure  6(a)).  The  biased  autocovariance  estimate  from  Eqn.  (5)  is  seen 

in  Figure  6(b)  to  be  piecewise  linear.  Natural  edges  tend  to  be  less  abrupt, 

resulting  in  smoother  autocovariance  functions  which,  over  small  lag  values, 

fit  the  general  shape  of  the  exponential  in  Eqn.  (1).  Thus,  Eqn.  (9)  is  seen 

to  behave  well  if  N  is  small  enough  to  traverse  single  edge-like  features. 

Edge  height  (contrast)  is  irrelevant  to  the  estimator.  Figure  7  shows  maps 

of  /  p  and  /  p  for  the  well-known  "girl"  image  for  the  cases  N  •  16,  and  32. 
1  2 

The  case  where  N  ■  16  is' near  optimum  in  terms  of  resolution.  (The  case 
H  *  8  is  excessively  noisy.) 


Block  estimation  of  parameters. 

In  certain  applications,  such  as  image  data  compression,  single  values 
of  the  autocovariance  parameters  «uce  sufficient  for  representing  blocks  of 
N  x  N  pixels,  rather  than  values  at  every  pixel.  We  may  think  of  the  esti¬ 
mation  procedure  as  before  with  an  additional  parameter  I  -  the  increment 
between  estimates  in  the  (i,j)  image  space.  (I  ■  1  in  all  the  preceding 
examples.)  The  estimation  of  /  p  values  is  then  posed  as  follows: 

•  Apply  Eqn s.  (4) ,  (8) ,  (9)  and  (10)  to  give  maps  consisting 


of  (D/I)  x  (D/I)  values  of  /  p  and  /  p  .  (D  is  the  image 

l  2 

dimension. ) 

•  Take  averages  over  blocks  of  dimension  (N/I)  x  (N/I)  to 

give  a  total  of  (D/N)  x  (D/N)  values  of  /  p  and  /  p  ; 

1  2 

i.e. ,  single  values  of  /  p"  and  /  p  for  each  N  x  N  -  pixel 
-  1  2 

subblock. 

Values  of  I  »  4  and  N  *  16  have  been  determined  by  trial  and  error  to 
/'.eld  results  which  correlate  well  with  observed  signal  activity.  Greater 


nr  resents  yield  insufficient  data  for  the  averaging  process. 


are 


Estima£ion_from_de2raded_Jma2es. 

The  need  to  estimate  autocovariance  parameters  from  images  which 
blurred  and  noisy  arises  in  the  process  of  image  restoration.  Experiments 
with  images  blurred  by  3  x  3  -  neighborhood  averaging  and  additive  Gaussian 
noise  (standard  deviation  S.O)  reveal  no  appreciable  change  in  the  maps  of 
z  or  V  p  ".  Even  noise  levels  of  two  or  three  times  this  produce  useable 
estimates  of  the  autocovariance  parameters. 


V.  DATA  COMPRESSION  APPLICATION. 

In  this  section ,  we  show  an  application  of  local  correlation  parameters 
to  block  interpolative  coding  £is3*  We  address  so-called  destination  intarpola 
tion,  where  the  rate  of  pixels  transmitted  from  image  lines  is  proportional 
to  /  p  .  The  procedure  is  as  follows* 


•  divide  the  image  into  N  x  N  -  pixel  subblocks. 


•  estimate  /  p^  for  each  subblock  (see  previous  section) . 


subsample  each  subblock  by  a  function  r  ■  X/V  p 

to  yield  a  subblock  of  N  x  N'  pixels,  (x  -  max  rail  /  p  -j  , 

1 

M*  -  N/r). 

transmit  the  subsampled  data,  plus  the  values  of  /  p  . 
reconstruct  the  image  by  interpolation  (the  inverse  of  sub¬ 
sampling. 


Subsampling  along  lines  is  convenient  in  images  which  are  raster  scanned. 
Additional  compression  is  possible  if  across-line  redundancy  is  reduced  by 


subsampling  in  this  direction  using  /  o  values. 

2 

Figure  8  shows  the  girl  image  after  compression  (8(a)  )  and  reconstruction 
(  8(b)  )  using  J  only.  The  compressed  image  consists  of  contiguous  blocks 
whose  width  is  proportional  to  J  .  Data  compression  is  only  30  %  ,  indica¬ 
ting  that  this  particular  data  contains  little  spatial  redundancy.  The  degra¬ 
dation  of  high  frequency  details  associated  with  nonadaptive  interpolative 
coding  is  considerably  reduced.  Incorporating  DPQ1  coding  would  reduce  the 
data  rate  by  an  additional  factor  of  two. 


VI.  CONCLUSIONS. 

A  simple,  useful  method  for  estimating  local  correlation  parameters  in 
Images  has  been  presented.  The  joint  assumptions  of  stationarity  within 
NxN-pixel  subblocks  (N»16)  and  the  nonseparable  exponential  model  for  auto- 
covariance  is  seen  to  work  well  when  judged  on  the  basis  of  comparisons  with 
observed  signal  activity.  The  resulting  maps  of  autocovariance  parameters  are 
useful  in  image  data  compression  and,  we  expect,  in  adaptive  filtering  and 


restoration. 


\  *. 
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